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Abstract: In this paper, a very simple method is used to derive the weakly 

singular traction boundary integral equation based on the integral relationships for 
displacement gradients [Okada, Rajiyah and Atluri, 1989J. The concept of the MLPG 
method is employed to solve the integral equations, especially those arising in solid 
mechanics. A moving Least Squares (MLS) inteipolation is selected to approximate the 
trial functions in this paper. Five boundaiy integral solution methods are introduced: 
direct solution method; displacement boundary-value problem; traction boundary-value 
problem; mixed boundary-value problem; and boundary variational principle. Based on 
the local weak form of the BIE, four different nodal-based local test functions are 
selected, leading to four different MLPG methods for each BIE solution method. These 
methods combine the advantages of the MLPG method and the boundary element 
method. 
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1. Introduction 


Numerical methods based on integral equation formulations of continuum mechanics, 
have become increasingly popular, for the solution of practical engineering problems. 
The integral equation methods can be considered to be derivable from the global 
unsymmetrical weak forms and Petrov-Galerkin approximation schemes (Zhang and 
Atluri, 1986, 1988). In these methods, the trial and test function-spaces are quite different 
from each other. The test functions correspond to fundamental solutions, in infinite space, 
of the differential operator of the problem, and hence are usually infinitely differentiable, 
except possibly at the singular points. If the fundamental solution can be derived tor the 
entire differential operator of the problem, the integral representations involve only 
boundary-integrals, whose descrctizations lead to the “boundary element method . 

In practice, information about the tractions on an oriented surface in the 
continuum is often required, so it is necessary to introduce an integral relation for the 
tractions. Such a traction integral equation is most often derived directly from the 
displacement integral equation, by a direct differentiation of the displacement integral 
equation. The resulting relation, however, contains a hypersingidar kernel. The numerical 
solution of the hypersingular traction integral equation is challenging; and various 
strategies have been developed to cojie with these difficulties. On the other hand, a 
distinctly different approach was presented by Okada, Rajiyah and Atluri (1989), 
[wherein they use the gradients of the fundamental solution as the test functions, instead 
of the fundamental solution itself], in order to obtain a singularity-reduced traction 
integral equation. Li and Mcar (1998), Li, Mear and Xiao (1998) also developed a 
weakly singular global wcak-form traction integral equation. They follow a systematic 
approach of regularization, in which stress functions associated with certain fundamental 
solutions are utilized to affect a global “inlegration-by-parts”. Based on the weakly 
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singular global weak-form integral equations, the symmetric Galerkin boundary element 
method (SGBEM) [Bonnet, Maier and Polizzotto (1998)] can be developed. However, in 
this paper, based on the singularity-reduced traction integral formulation of Okada, 
Rajiyah and Atluri (1989), we will develop a weakly singular local weak-form traction 
integral equation, without the introduction of stress functions. Our derivation is much 
simpler and direct than that of Li and Mear (1998). 

Meshless methods, as alternative numerical approaches to eliminate the well- 
known drawbacks in the finite clement and boundary clement methods, have attracted 
much attention in recent decades, due to their flexibility, and due to their potential in 
negating the need for the human-labor intensive process of constructing geometiic 
meshes in a domain. Such mcshless methods are especially useful in problems with 
discontinuities or moving boundaries. The main objective of the meshless methods is to 
get rid of, or at least alleviate the difficulty of, meshing and remeshing the entire 
structure; by only adding or deleting nodes in the entire structure, instead. Meshless 
methods may also alleviate some other problems associated with the finite element 
method, such as locking, element distortion, and others. 

Recently, the mcshless local Pctrov-Galcrkin (MLPG) method, has been 
developed, in two of its alternate forms, in Zhu, Zhang and Atluri (1998a, b), Atluri & 
Zhu (1998a, b) and Atluri, Kim and Cho (1999), for solving linear and non-linear 
boundary problems. This method is truly mcshless, as no finite elcment/or boundary 
element meshes are required, cither for the purposes of interpolation of the trial and test 
functions for the solution variables, or for the purpose of integration of the ‘energy’. All 
pertinent integrals can be easily evaluated over over-lapping, regularly shaped, domains 
(in general, spheres in three-dimensional problems) and their boundaries. Remarkable 
successes of the MLPG method have been reported in solving the convection-diffusion 
problems [Lin & Atluri (2000)]; Iraclurc mechanics problems [Kim & Atluri (2000), 
Ching & Batra (2001)]; Navicr-Stokes flows [Lin & Atluri (2001)]; and plate bending 
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problems [Gu & Liu (2001); Long & Atluri (2002)]. A comparison study of the 
efficiency and accuracy of a variety ol mcshless trial and test functions is presented in 
Atluri and Shen (2002a, b), based on the general concept of the meshless local Petrov- 
Galerkin (MLPG) method. 

In summary, the MLPG method is a truly meshless method, which involves not 
only a meshlcss interpolation for the trial functions (such as MLS, PU, Shepard function 
or RBF), but also a mcshless integration of the weak-form (i.e. all integrations are always 
performed over regularly shaped sub-domains such as spheres, parallelopipeds, and 
ellipsoids in 3-D). In the conventional Galerkin method, the trial and test functions are 
chosen from the same function-space. In MLPG, the nodal trial and test functions can be 
different; the nodal trial function may correspond to any one of MLS, PU, Shepard 
function, or RBF types of interpolations; and the test function may be totally different, 
and may corre.spond to any one of MLS, PU, Shepard function, RBF, a Heaviside step 
function, a Dirac delta function, the Gaussian weight function of MLS, or any other 
convenient function. Furthermore, the physical sizes of the supports of the nodal trial and 
test functions may be different. These featuies make the MLPG method very flexible. 
The MLPG method, based on a local fonnulation, can include all the other meshless 
methods based on global formulation, as special cases. 

In this paper, we use the concept of the MLPG method to formulate 
computational approaches to solve the integral equations, especially those arising in solid 
mechanics. A Moving Least Squares (MLS) interpolation is used to approximal the trial 
functions in this paper, although this is arbitrary. Four different nodal-based local test 
functions are also selected, leading to four dilferent MLPG methods for BIE. Based on 
the MLPG concept [Atluri and Shen (2002a)j, we label these variants of the MLPG 
method for solving the boundary integral equations as MLPG/BIEI, MLPG/BIE 2, 
MLPG/BIE 5, and MLPG/BIE 6, respectively [in conformity with the labeling introduced 
in Atluri & Shen (2002a, b)]. 
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The paper is organized as follows. In Section 2, we introduce the MLS 
approximation over the surface of a 3-D solid, using arbitrary curvilinear co-ordinates. In 
Section 3 are derived the weakly singular boundary integral equations, based on the result 
from Okada, Rajiyah and Atluri (1989), and the results are shown to be identical to those 
in Li and Mear (1998). In Section 4, five boundary integral solution methods are 
introduced: the direct solution method; displacement boundary- value problem; traction 
boundary-value problem; mixed boundary-value problem; and the boundary variational 
principle. Local weak forms for the boundary integral equations are introduced in each 
of these boundary integral equation (BIL) solution mcW\od&-, four different nodal-based 
local test functions are selected; and [hus four different MLPG methods are developed in 
each of those BIE solution methods. The paper concludes in Section 5. 

2. Moving least-squares method (MLS) in Curvilinear Coordinates, on 
the Boundary of a 3-D Body 

Since the nodes lie only on the boundary dQ (a curved surface in 3-D space) of a 3-D 
body curvilinear co-ordinates are necessary to define the MLS interpolates on the 
surface. The moving least-square method is generally considered to be one of the best 
schemes to interpolate data with a reasonable accuracy. The MLS interpolation does not 
pass through the nodal data. Here we give a brief summary of the MLS approximation for 
curvilinear co-ordinates. For details of the MLS approximation, see Belytschko et al. 
(1996), Atluri, Cho and Kim (1999), and Atluri and Shen (2002a). 

Consider a surface Fx, which is defined as the neighborhood of a point s and 
denoted as the domain of definition ol the Ml^S approximation for the trial function at s, 
and which is located in the pioblem surface dii. To approximate the distribution of the 
function u in Fx, over a number of randomly located nodes {s/}, /=1, 2,..., N, the moving 
least squares approximant u'{s) of u, VseFx, can be defined by 
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«''(s) = p' (s)j(s) VseF* 


( 1 ) 


where p\s)=[/;>i(s), ;; 2 (s), Pm(s)] is a complete monomial basis, and a(s) is a vector 
containing coefficients cfy(s), /=!, 2, which are functions of the spatial cuiyilinear 

coordinates s=[,v', The commonly used bases in 2-D are the linear basis, due to their 
simplicity. The linear basis assures that the MLS approximation has the linear 
completeness. Thus, it can reproduce any smooth function and its first derivative with 
arbitrary accuracy, as the approximation is refined. It is also possible to use other 
functions in a basis. For example, in problems with singular solutions, singular junctions 
can be included in the basis. 

The coefficient vector a(s) is determined by minimizing a weighted discrete Lj 
norm, which can be defined as 


(s) = S M-'/ (s / “ ' f 

/ = ! 


( 2 ) 


where wy(s) is a weight function associated with the node /, with w/(s)>0 for all s in the 
support of vv’/(s), S/ denotes the value of s at node /, N is the number of nodes in Tx for 
which the weight functions w/{s)>0. 

Here it should be noted that ii' , /=!, 2,..., N, in equation (2), are the fictitious 
nodal values, and not the actual nodal values of the unknown trial function m'’(s). 

Substituting a(s), which is solved by minimizing J in equation (2), into equation 
(1), give a relation which may be written in the form of an interpolation function similar 
to that used in the FEM, as 

u'’{s)='^<p' {s}i' ■, u’'{s,) = It' ^ u' , se I’x (3) 

/=i 
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where 


,/=i 


(4) 


where 


A(s)=^iv/(s)p(s/)|’'(s/) 

/-I 


(5) 


B(s) = [lV| (s^j(s| ), U'2 (s)l)(s2 ), • • • , (s)p(s;^, )] 


( 6 ) 


The partial derivatives of <p'{s) are obtained as 

M 


where A“' = (A“' denotes the derivative of the inverse of A with respect to a*, which is 
given by 


A;=-A“'A,A' (8) 

The MLS approximation is well defined, only when the matrix in eqn (5) is non- 
singular. 0^(s) is usually called the shape function of the MLS approximation, 
conesponding to the nodal point x/. From eqns (4) and (6), it may be seen that 0(s)=O 
when w/(s)=0, that preserves the local character of the moving least squares 
approximation. The nodal shape (unction is complete up to the order of the basis. The 
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smoothness of the nodal shape function is determined by that of the basis, and of the 
weight function. 

The choice of the weight function is more or less arbitrary, as long as the weight 
function is positive and continuous. Both Gaussian and spline weight functions with 
compact supports can be considered in the present work. The Gaussian weight function 
corresponding to node 7 may be written as 


w,{s) = 



Q<d, < r, 
d, >r, 


(9) 


where t/y =|s-s/| is the distance from node sy to point x, cy is a constant controlling the 
shape of the weight function vvy (and therefore the relative weights), and ry is the size ol 
the support for the weight function U 7 (and thus determines the support of node S/). 

A spline weight function is defined as 


w 


/(x) = 


1-6 


^/,V 


r y 




V ' y 


0 , 


V y 






V ' y 


, 0<d, <r, 

d, >r, 


( 10 ) 


The size of support, ri of the weight function rvy associated with node 7 should be 
chosen such that /■/ should be large enough to have a sufficient number of nodes covered 
in the domain of definition of every sample point {ii>in}, in order to ensure the regularity 
of A. It can be easily seen that the spline weight function (10) possesses C continuity. 
So, the MLS shape functions, and the corresponding trial functions are C' continuous 
over the entire domain. 
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3. Boundary Integral Equations in Linear Solid Mechanics 

Let Oij be the Cartesian components of Ihc Cauchy stress tensor, and let fi be the body 
force per unit volume. The equations o( linear and angular momentum balance are: 

ay,.,+/, =0; (11) 

where ( )„ denotes differentiation with respect to material coordinates .r,-. For a linear 
elastic isotropic solid, the stress-strain relations are: 

^ij == ( 1 ^) 




where, A and jU are Lame constants and Sij is the Kronecker’s delta. The strain- 
displacement relations aie 


^ki ' ''''/.*■) 


( 14 ) 


The boundary tractions arc given by: 


where /j, are components of a unit outward normal to the closed boundary dQ=F. 

Let Ui be the trial functions for displacement, and let u. be the coiresponding test 
functions. The global weak-form of the equilibrium equation (1 1), can be written as: 
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Lk , +/>/'“ =0 


(16) 


where Oij are assumed lo be wriUen as functions of m, through eqns (12) and (13). 

We assume that the test functions «, are the fundamental solutions, in infinite 
space, to the Navier equations of elasticity, i.e. 

+f5(x-^ =0 (17) 

where Cp denotes the direction ol the unit load at We assume that the solid is 

isotropic, in which case the solution (or (17) is readily available. Thus, 

(nosum for /;) (18) 


and 


r, 09) 

Here, is the yth component of displacement at location x,„ due to a unit load along the 
plh direction at location Likewise, is the yth component of traction on an oriented 
surface at x,„ due to a unit load along the pth direction at By using the divergence 
theorem twice on eqn (16), and substituting eqn (17) in the resulting equation, one 
obtains the integral equation: 

where T is the global boundary; and is the global domain, enclosed by F . By taking the 
point in the limit, to the boundary, one may obtain the well-known boundary integral 
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equation for u,,. The kernel f], is singular at the source point with the order of 
singularity being 0{l/r^) in 3D or 0{\/r) in 2D; while the kernel w*^ is weakly singular, 
with the order 0{\/r) in 3D or 0(ln(J//)) in 2D. Hence, the order of the singularity of cqn 
(20) is D(l/r^) in 3D or 0(l/r) in 2D. 

It is also well known that the singular kernels uj^, and remain integrable in the 
limit when tends to the boundary. By a direct differentiation of eqn (20) with respect 
to one may obtain; 


di( 






duy{\.^) 




,/ :(x)- 


c>yx,^) 




L/r+f — -(in 

1 JnDV f 


(21) 


In the limit as ^ -^F, the kernel is singular with the order of singularity being 

D(l/r^) in 3D; the kernel becomes hyper-singular, with the order of singularity 

being 0(l/r^) in 3D, and thus becomes numerically intractable. Hence, eqn (21) becomes 
hyper-singular with the order of singularity being 0(l/r^) in 3D. 

To circumvent these difficulties, Okada, Rajiyah and Atluri (1989, 1994) 
developed alternate types of integral rcpiesentation for displacement gradients, which 
will be introduced in the following, while they originally presented the non-hypersingular 
integral equations for the Field-Boundary Element Method in nonlinear solid mechanics. 

Instead of writing the global weak-form of the linear momentum balance relations 
in a scalar-form as in eqn (16), Okada and Atluri (1989, 1994) wrote the global weak- 
foiTns of the linear momentum balance relation in a three component vector form, as: 

Lk ,„+// K .''«=0 

Assuming that the linear clastic solid is homogeneous, (22) can be rewritten as: 
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(23) 



+ 


/. ]uj,dQ = 0 


Integrating eqn (23) by parts, applying the divergence theorem, and making use of eqn 
(15), one obtains: 







m m,k 


E. ii )f/r+f Ju,.dQ.+ \ f.S.,dQ = 0 

: ijfnn m,n j,i ^ Ju ^ ^ ^ 


(24) 


where 


f, 





), 


(25) 


and 




(26) 


Upon taking «. and !■ to be the fundamental solutions as in eqns (18) and (19), one 
obtains the global integral relation: 

= J,, ('b‘< ,<^Lp - 'h,><.„.k^Lp - t/jp,k w - f/jp.k^^ 

In eqn (27), C=1 when and C=l/2 at a smooth part of the boundary. As compared 

to eqn (21) wherein the kernel is involved in the boundary integral, in the 

presently derived eqn (27), only the kernels // and are involved at F. Note that the 
orders of singularity in w*, . and , which aie equal, are nevertheless smaller than that 
in dt*ji,jd^^ . The singularities in the integrals in eqn (27) are in general tractable, and the 
integrals in eqn (27) can be evaluated by the method suggested by Guiggiani and Casalini 
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(1987) for two-dimensional case, by Guiggiani and Gigante (1990) for three-dimensional 
case, or by an alternate method (Okada and Atluri, 1994). Once the tractions tj at the 
boundary are completely known, the displacement gradients Uij can be determined from 
eqn (27) in the interior Q, as well as at f. 

Now, we will further regularize ecjn (27), in order to reduce its singularity to be 
same as the orders of singularity in »* , . The surface curl operator can be defined as 


where is the usual permutation symbol. The surface curl, eqn (28) is associated with 
the following form of the Stokes’ formula: 

(f).,/(x)/r(.)=£/(x>(x. (29) 

where F is any regular surface with the edge contour 3F. 

The first two terms in the extreme right side of eqn (27) can hence be written as 

nmp ~~ mnp ~ ^ nwp ^ nmp^ nk( ^ I m ( 30 ) 


Thus, eqn (27) can be rewritten as 


a 




Jii ( 31 ^ 

)r/r - [|, fjuln.dT - £ /. ,«;rfn(x)] 


The orders of singularity in and cr,*,„y, , which are equal, are 0(l/r ) in 3D. So, eqn 
(31) is singular with the order of singular G(l/r^) in 3D. This equation (31) is exactly the 


14 



same as that in Li, Mear ami Xiao (1998). However, here we used a totally different and 
simple method to derive it, based on the traction integral equation developed by Okada, 
Rajiyah and Atiuri (1989). In Li, Mear and Xiao (1998), a multi-valued stress function is 
introduced, and the derivation is based on the gradient of the displacement integral 
equations. 

An integral relation for the traction vector then follows immediately from eqn 
(31) by an application of Hook’s law, with the result: 

-[XP,i,u,P99n ( 32 ) 

where -^) = . It is noted that the symmetry condition 

(x -<^)= mV (x - 1 ^) has been used. Eqn (32) is singular with the order of singular 
0(l/r^)in 3D. 

3.1 Weakly singular integral equations 

Using Somigliana’s identity, a weakly singular displacement integral equation can 
be readily obtained. This can be achieved by a subtraction technique often motivated by 
consideration of a rigid body motion (Brebbia and Dominguez, 1992). Li, Mear and Xiao 
(1998) exploited the decomposition of the stress fundamental solution, in order to obtain 
a weakly singular displacement integral equation. Here, we will also use this method. The 
stress fundamental solution can be decomposed as 

(x - «) = <’,..,0',:; . (>1 - P- 4 ) ( 33 ) 


where 
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(34) 


M (l - v)G ,1' (z) = (l - ^ ^ 


and 




z. 


(35) 


where v is the Poisson’s ratio and r(z) = ||z| . 

Similar to Li, Mear and Xiao (1998), the following results will be used: 






Irt 


a.v. 


+ 


dx 


(36) 


where 




2^^S,„ + 2(1 - + 4v5,.5„ - 2<\,S 


(37) 


and 


8^(l-^)/u,(z) = y 


(l-2rK - 


^kji + 


(38) 


where fj, is the shear modulus. Thus, the traction global integral equation can be obtained 
from eqa. (35-38) as [same as in Li, Mcar and Xiao (1998)]: 
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Ct, {^)=n, ^ G/„, (x - 0 , (x >/l ' + n „ ~ J_. C (x - ^ )Ojt j (x>/r 






( 39 ) 


Incorporating the decomposition (21) into (20), and utilizing Stokes’s tlieorem to 
integrate by parts, a weakly singular integral equation is then obtained as 


(d-; (x.«)'r + ^ J/., 


+ 


I, Dji, (x,^>/r + £ fj (x>; {x,^)dn 


(40) 


Note that the kernel /i, (x)(<^, - a, )//-^ appearing in the second integral is only weakly 
singular since «, (x)(^,- -.v, )//•-> 0(/) as r^O, hence, eqn (40) involves only 6>(l/r) 
singular kernel. 

On the other hand, we can find that 


^k^ipnn^m.n^ ip.i ^ jl’.k 

= ('h"]p., - >V‘l,k ) = (^.J^in,k ^^n^Jp 


(41) 


Thus, using eqn (41), the traction integral equation (39) can also be rewritten as 


Cl, {4) 


~ ^ithpk J|. (^ij^'imk jp Knp^^rn.k 

- ^nl,pk I j|, fj » I'h - £ fj.k li]p 


(42) 


The corresponding weakly singular integral equation can then be obtained as 
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( 43 ) 



Note that the singular kernel appearing in the integral equation (43) is only weakly 
singular (Cauchy singular kernel), with the order of singular C(l/r), similar to the 
displacement integral equation (40) and weaker than that in equation (39). In this 
equation, second derivatives of the shape 1 unctions arc needed in constructing the global 
stiffiness matrix. It is important to note that even though the singularity of the integral 
kernels has been reduced to Cauchy type, it still requires that the derivatives of the 
displacements are continuous. However, the calculation of the derivatives of shape 
functions from the MLS approximation is quite costly. Hence, we choose to use equation 
(39) to develop the MLPG approach in this paper. 

4. Solution Methods for BIE: the Meshless Local Petrov-Gaierkin 
(MLPG) Approach 

To solve the boundary integral equations, different methods are proposed in this section, 
similar to those in traditional boundary clement methods (Atluri and Grannell, 1978). 

To start with, we note that in a well-posed boundary-value problem in 3-D solid 
mechanics, at any point on the boundary, any 3 quantities out of the 6, viz, w, (three) and 
tj (3) can be prescribed and the other 3 components are the unknowns to be solved for. In 
a pure displacement boundary-value problem, all the 3-components of Ui are prescribed at 
every point on the surface of the 3-D solid; and correspondingly, all the three components 
of tj at every point on dQ are the unknowns. Conversely, in a pure traction boundary- 


18 




value problem, all the three ij are prescribed at every point on dQ, while all the 3 
components of iii at 3Q arc unknown. 


4,1 Direct solution method 

This method attempts to obtain a direct numerical solution to equation (40). Upon 
multiplying (40) by a continuous test function (jp(^) over a local sub-domain Ts (a part of 
T) and applying the divergence theorem, a local ''weakly singular' weak-form of the 
global displacement integral equation lor a system, without body force, is obtained as 

+ :rJr (44) 

To obtain the discrete equations based on MLS intejpolation, the following 
interpolations are used 




/-) 


(45a) 


(45b) 


^y,(s)=S'/^'(s>y,' 


/-I 


(45c) 
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where s is a point with curvilinear co-ordinates (^ 1 ,^ 2 ). 0^(s) and y/'{s) are the nodal 
shape functions for trial and test functions centered at node, ,V| and sj are functions of .v,. 
In general, in meshless interpolations, ul are fictitious nodal values. Substitution of eqns 
(45) into formulations (44) leads to the following discretized system of linear equations 
for each node / 

m 

where 

^lp.o = 1,;' V/' 0' (x><;(x.0^r4r, (47) 

u T„jj = V/ ' J, ( // (x - + g,;; (x - ^ )d„^ ^ (48) 

w„ c<j^//'4r, (49) 

m is the number of the background mesh, and / is the nodal number. Due to the fact that 
the first step of the dual integration is performed on the global boundary surface F, a 
background mesh will have to be used on F. It is noted that the repeated indices imply 
summation here. The trial functions appear only in the global integration over F, which 
means that the integrand for the local integration over Fs will be simple. It can be found 
that the assembly process is not required to form a global ‘stiffness’ matrix. 

The equations can be solved in the same way as in the conventional BEM, except 
that the transformations between and /T/ , f/ and // must be performed, due to the 
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fact that the MLS inlcipolales lack the della function property of the usual BEM shape 
functions [Atluri, Kim and Cho (1999), Atluri and Shen (2002a)]. 

On the boundary surface on which m, arc prescribed, h/ can be obtained by tlie 
following transformation: 

u’ =^K,jn^ (50) 

J=\ 


On the boundary surface on which /, are prescribed, // can be obtained by the same 
transformation as the above 

(51) 

, 7-1 

where R,j =i^'^(s/)] '■ The details of the transformation method for the imposition of 
boundary conditions can be lound in [Atluri, Kim and Cho (1999), Atluri and Shen 
(2002a)]. 

Based on the local weak form, the meshless local Petrov-Galerkin (MLPG) 
method is used to treat the global boundary integral equations. However, due to the fact 
that a global weak form is used to derive the boundary integral equations, we need a 
background mesh to perform the global integration first. For the second step, to perform 
the local integration, no mesh is used. 

In general, in MLPG, the nodal trial and test functions can be different, the nodal 
trial function may correspond to any one ol: 1. MLS; 2. Partition of unity; 3. Shepard 
function; or 4. Radial basis function interpolations (Atluri and Shen, 2002a); and the test 
function may be totally different. Furthermore, the size of the sub-domains over which 
the nodal trial and test functions are, respectively, non-zero, may be different. Different 
choices for the basis functions for the trial function, and the test functions, will lead to 
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different approximation methods. We will use the MLS interpolations introduced in 
Section 2, to generate the trial lunctions. Ba.scd on the concept of the MLPG, the test 
functions over each local boundary can be chosen through a variety of ways. 

As a known test function is used in the local weak form (LWF), the use of the 
LWF for one point (and here for one sub-surface F,) will yield only one algebraic 
equation. It is noted that the trial functions u or t within the sub-surface F,, in the 
interpolations without Kronecker’s Delta properties, is determined by the fictitious nodal 
values u' or t' , respectively, within the domain of definition for all points x falling 
within F.v. One can obtain as many equations as the number of nodes. Hence, we need as 
many local boundary surfaces F, as the number of nodes in the global boundary surface, 
in order to obtain as many equations as the number of unknowns. In the present paper, 
dealing with 3D problems, the local boundary surface is chosen as a circle, centered at a 
node X/. Four test functions are chosen to lormulate four different MLPG methods for 
BIE: 

(1) the test function over F, is the same as the weight function in the MLS 

approximation: The resultant Meshless Local Petrov-Galerkin Method for directly 
solving the BIE is denoted as MLPG/BIEI. In this case, just let in eqns 

(47) , (48) and (49). 

(2) the test function over F., is the collocation Dirac’s Delta function (collocation 
method): The resultant Meshless Local Petrov-Galerkin Method for directly 
solving the BIE is denoted as MLPG/BIE2. In this case, let y/=5{^') in eqns (47), 

(48) and (49), then it can be Ibund that we get the collocation BEM. The 

collocation BEM can be treated simply as a special case of the MLPG approach. 
Chati, Mukheijce and Mukheijec (1999) developed a boundary node method by a 
coupling between boundary integrals equations (BIE) and Moving Least-Squares 
(MLS) inteipolants. In fact, their method can be treated simply as a special case of 
the present MLPG method approach lor BIE, if we let in the local weak- 


22 


form. However, here we use cm entirely different regularization approach, to 
avoid the strong singidarity. 

(3) the test function over T, is the Heaviside step function (constant over each local 
sub-domain Fv): The resultant Mcshless Local Petrov-Galerkin Method for 
solving the BIE is denoted as MLPG/BIE5. In this case, let i//=l in eqns (47), 
(48) and (49). 

(4) the test function over F, is identical to the trial function (Galerkin method): The 
resultant Meshless Local Petrov-Galerkin Method for solving the BIE is denoted 
as MLPG/BIE6. In the Galerkin method, the trial and test functions come from 
the same space. In this case, let r/, =»,(x), i.e. let i//=^ in eqns (47), (48) and (49). 
The interrelationships of these developments can be illustrated as in Fig. 1. 

Underlying all these mcshless methods for BIE is the general concept of the mcshless 
local Petrov-Galerkin method; thus, MLPG provides a rational basis for constructing 
meshless methods with a greater degree ol flexibility. Theoretically, as long as the union 
of all local domains covers the global boundary surface, the boundary integral equation 
will be satisfied. 

In the collocation BEM, the location of collocation nodes is an important 
ingredient for the success of the method. However, if we employ MLPG method, this 
issue can be avoided. 

4.2 Displacement Boundary- Value Problems in Solid Mechanics 

In this case, displacements are prescribed on all the boundaries. Then, the local weakly 
singular weak-form displacement integicd eciuation (44) for a system without body force 
can be rewritten as 
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( 52 ) 


j,. <7/, )j|. I j{^)* ]p ^ V = j| f }lp (^ V'^T 


To obtain the discrete equations based on MLS interpolation, using interpolations (45b, c), 
we can obtain 



}l) 


(53) 


where all (j are unknowns, and 

( 54 ) 

//4(x,0'„,(0L(x)/r4r, (55) 

/L,'T/(x)r;;.(x-0/r4r, 

Similar to Subsection 4.1, four test functions are chosen to formulate four different 
MLPG methods for the displacement 13 IE: 

(1) MLPG/DBIEl: the test function over E, is the same as the weight function in the 
MLS approximation. In this case, just let i//=w^in eqns (53)-(55). 

(2) MLPG/DBIE2; the lest function over E., is the collocation Dirac’s Della function 
(collocation method). In this case, let i//=<5(^^) in eqns (53-55), then it can be 
found that we gel the collocation displacement BEM. 
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(3) MLPG/DBIE5: ihc test function over T, is the Heaviside step function (constant 
over each local sub-domain r,v). In this case, let i//=l in eqns (53)-(55). 

(4) MLPG/DBIE6: the test function over I',, is identical to the trial function (Galerkin 

method). In this case, let t/, =h(x)> i-^- i'l eqns (53)-(.‘>5). 

The interrelationships of these developments can also be illustrated as in Fig. 1. 


4.3 Traction Boundary-Value Problem in Solid Mechanics 


In this case, tractions arc prescribed on all the boundaries. Upon multiplying (39) by a 
continuous test function viXO and applying the divergence theorem, a local weakly 
singular weak-form traction integral equation is obtained as 


= -fer ' I,. ^r'',(0jpG,,„y(x-^)D,„H^(x>/r4r, 


(56) 


where dFj is the edge of the local boundary surface Fj. It is noted that only the weakly 
singular kernel uj^, is involved in these equations. In the present formulations, the 
boundary integral equations are satislied in all the local boundary surfaces Fj. 
Theoretically, as long as the union of all local boundary surface covers the global 
boundary surface, i.e., LJEgDr, the boundary integral equations will be satisfied in the 
global boundary surface. 

To obtain the discrete equations based on MLS interpolation, the following 
inteipolations are used 

",(s)=if(pi/ (57a) 
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(57b) 


/-I 

Substitution of eqns (57) into formulations (56) leads to following discretized system of 
linear equations for each node / 

(58) 

m 


with all Uj are unknowns, and 


H 


ni 

IpJj 


- £,./ V ' j| (x -e)Dj ■' {xyivd^, 


(59) 


2 /;, =-J,,cf;,(^V^/r; o„w'l G/„,{’^-^rj{x)drdi\ 

+ V ' J,- g4> Jj ( 60 ) 


Similar to Subsection 4.1, four test functions are chosen to formulate four different 
MLPG methods for the traction BIE: 

(1) MLPG/TBIEI: the test function over r, is the same as the weight function in the 
MLS approximation. In this case, just let y/=w' in eqns (59) and (60). 

(2) MLPG/TBIE2: the test function over r, is the collocation Dirac’s Delta function 
(collocation method). In this case, let \j/=5{^') in eqns (59), and (60), then it can 
be found that we get the collocation displacement BEM. 

(3) MLPG/TBIE5: the lest function over Py is the Heaviside step function (constant 
over each local sub-domain Py). In this case, let i//=l in eqns (59), and (60). 
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(4) MLPG/TBIE6: the test function over r, is identical to the trial function (Galerkin 
method). In this case, let v, =Ui(x), i.e. let in eqns (59), and (60). 

The interrelationships of these developments can also be illustrated as in Fig. 1. 

4.4 Mixed Boundary-Value Problem in Solid Mechanics 

This method is developed for the special problems of fracture mechanics of linearly 
elastic three-dimensional solids, containing cracks. The formulation is based both on the 
displacement integral ecjuation (40) and traction integral equation (39). When applied to 
fracture problems, the traditional boundary clement methods (or the direct methods in 
section 4.1) become mathematically degenerate in that when the displacement integral 
equation is applied to points on the crack surface, information about the traction on the 
crack is lost (Cruse, 1988). To circumvent this difficulty an integral equation for the 
traction on the crack surface may be employed. Moreover, in practice, information about 
the tractions on an oriented surface in the continuum is often required, so it is necessary 
to use both the traction and displacement inlcgral equations. 

Upon multiplying (40) by a continuous test function c]p{^ and applying the 
divergence theorem, a local weakly singular weak-form of the displacement integral 
equation for a system without body force is obtained [as in (44)] as 

|q-, te>/,fe)rfr. = 1,^ , (x v; (x.^ >/r<r, 

+ (61) 

•'I r 

+ 1 , 
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Upon multiplying (39) by a continuous test function Vb(^) and applying the 
divergence theorem, a local weakly singular wcak-form of the traction integral equation 
is obtained as 


-jQ,fe>X«Vr, =1 ov.(x-^),(x>rrfr. 


■£/, («)}, c,, (x - (xVrrf^, - v,g)..e)[ 


where dFs is the edge of the local boundary surface Fs. It is noted that only the weakly 
singular kernel is involved in these equations. 

In the absence of body forces, let the regular boundary be partitioned into a 
portion F, on which tractions are preseiibed and a portion F„ on displacements are 
prescribed such that F=F, +F„. Applying the local weak-form displacement integral 
equation (61) on F„ with U, and, similarly, applying the local weak-form traction 

integral equation (62) on F, with i’r=0 on F„ gives rise to the formulation 

A{q,t)+ B{q,u) = F{(j) (for F„) (63) 


G{v,t)+Uiv,i{} = Q{v) (foiF,) 


(64) 


where 


( 63 ) 
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‘jpj/<x) ''h'V/( 5 - x) '"-'j j( P) ' 


C- 



W*) ’■) “w J(?) "* 1 



where « and t are Ihe prescribed dis|)lacement and tractions, respectively, on the 
boundary Fu, and on tbe boundary T,. 

To obtain the discrete equations based on MLS interpolation, using the 
interpolations (45) and (57), and substituting them into formulations (63) and (64) leads 
to the following discretized system of linear equations for each node 1 

Sw:v7 +««»;)='■» CD 

m 

E(cv7+";;vip=e,, ov 

m 

where 








(75) 


dVcir 


H 


m 

fpJj 




(76) 
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/■V =I,i -j,, (*)<;(x.'5)'rar. 

A.i<,(x)c,;;(x-|);rar, 

* Sll II 


( 77 ) 


+ £,, v ' £ o,; (x - ? X, (li >'i '‘(s. 

+|,,v''".te)|,, -^);(4'ivr, (78) 

-£,d,v,'£ c„(x^^yj.,T,(x>/r</r, 

+ £, w'l_ c„„ (x - 5 Y\u, (x)7r () 5 , 


m is the number of the background mesh, and / is the nodal number. It is noted that the 
repeated indices imply summation here. Note again that the trial functions appear only in 
the global integration, which means that the integrand for the local integration will be 
simple. 

Similar to Subsection 4.1, four test functions are chosen to formulate four 
different MLPG method.s for the mixed BIB: 

(1) MLPG/MBIEl; the test function over B, is the same as the weight function in the 
MLS approximation. In this case, let \j/=vv in eqns (73), (74) and (77), and 
rewrite eqns (75), (76) and (78) as 


-I G;(x-^>V/rr/i; (x-^>Vrr/r, (79) 




( 80 ) 
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( 81 ) 


<2„ = cf„ )vv, cir^ - D„ w, g; (x - ^ )■ . {x)dr dr, 

- j,., ^ 'G I, (x - ^ ; (x ; 

Here that the test function (weight function) vanishes at is considered. 

(2) MLPG/MBIE2: the test function over Fv is the collocation Dirac’s Delta function 

(collocation method). In this case, let in eqns (73), (74) and (77), then it 

can be found that we get the collocation BEM. 

(3) MLPG/MBIE5: the test function over F,, is the Heaviside step function (constant 
over each local sub-domain F.,). In this case, let i//=l in eqns (73), (74) and (77), 
and eqns (7.5), (76) and (78) can be rewritten as 




(82) 


// 


m 


|,,„C-,„„(x-^)D„0^(x>/Fr/^, 


(83) 


Q„,=-[,ct„{^)dr, +£, £6;;(x-^]F;(xK^^™ 

+ [, ^h,(^)|,./2,;;(x-0“ (x>/rr/F, (84) 

+ j, C,,,,,, (x - >7 ,,7, (x>/F d^, 


It is well-known that the numerical integration plays an important role in the 
convergence of numerical solutions of mcslilcss methods. From equation (83), it can be 
seen that some surface integrals over F,, replaced by the curve integrals over dFs, which 
will improve effectiveness of this method. 
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(4) MLPG/MBIE6: the lest funclion over f, is identical to the trial function (Galerkin 
method). In the Galerkin method, the trial function and test function come from 
the same space. In this case, let <^,=;,(x) and v,=Mj(x), then we have 

(85) 


= I,, f - «>-„)* 'rfr</r, + )<'r,/r, (86) 


n„ =i; QT„(s'>'jr -l,yli,(xy^^{x.(}irdr, 

(89) 


e„ =-|^,o,g>'</r. G;(x-^)f,(x>rrfr, 

+ [,f„.g)(//i;(x-^;r,(xVr</r, (9o> 

- j,; I, c,,,,(x-|>VT,(x>/rar, 
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If we replace the local boundary by the global boundary in the above equations, the 
foiTnulations are same as those in SGBEM [Li, Mcar and Xiao (1998)], 

The interrelationships of these developments can also be illustrated as in Fig. 1. It 
is noted that in all the MLPG methods in this study, the usual “element assembly’’ 
process is not required, unlike the SGBEM, to form the global stiffness matrix. 

Due to their flexibility, and due to their potential in negating the need for the 
human-labor process of constructing meshes along boundary surfaces, such MLPG 
methods for BIE are especially useful in those problems with discontinuities or moving 
boundaries. The main objective of the mcshless methods is to get rid of, or at least 
alleviate the difficulty of, meshing and remeshing the entire boundary surface; by only 
adding or deleting nodes in the entire boundary surface, instead. 

The MLPG methods lor BIE are characterized by weakly singular kernels and 
meshless local weak form. SGBEM is also characterized by weakly singular kernels, but 
it is based on the global weak form and meshes [i.e., in SGBEM, two global integrals are 
involved]. 

In SGBEM, for a specific pair of elements, which have no common points, the 
ordinary Gauss rule can provide sufficient accuracy. It is not appropriate for other 
elements, which are coincident or have one edge or one vertex in common. To deal with 
weakly singular integrals, there exist methods using transformation of variables in order 
to weaken or cancel out the singularity by Jacobian of the transformation before applying 
the ordinary Gauss rule. Thus, in SGBEM, there exist 9 ‘weakly singular element , 
which are coincident or have one edge or one vertex in common with the master element, 
for every master element for quadrilateral boundary elements, and 13 ‘weakly singular 
elements’ for every master element for triangular boundary elements, that need the 
transformation of variable. This procedure will be costly. However, in the MLPG for 
BIE, for every (master) local boundary surface Es, we can choose an appropriate 
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background mesh and adjust tlie size of Ps, to make Ps only involved in one background 
element, that will improve the cITectiveness obviously. 

In fact, this mctliod can be derived directly from the local weak-form of the 
boundary conditions: 


j,_ (», --zk</r=o 


(91) 


I, (', -(■, >,</r = o 


(92) 


where, a priori, the equilibrium equations and constitutive relations are satisfied. Using 
equations (40) and (39) to repiescnts m, and /, respectively, we can obtain the method 
developed in this subsection. 

4.5 Boundary Variational Principle 

This subsection deals with boundary solution methods, based on a boundary variational 
principle [Atluri and Granncll (1978)]. The functions that are assumed to be independent 
are: displacement u in the domain; boundary displacement u ; and boundary traction 
t [which is a Lagrange multiplcr to enforce u=u at P]. 

The corresponding vaiiational functional Pip, for linear elasticity, is defined as 
follows: 
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where the boundary displacement u satisfies the essential boundary condition, i.e., u = u 
on Fu, is the prescribed traction on the traction boundary F,, and no body force is 
considered. 

By carrying out the variations it can be shown that 


■5n, = f 


o 


V.,/ 


Sll, clQ + 


I (/, - 1. )5u.cir - {u. - U. )St_dr (94) 


with the vanishing of 6ri|,, one can also have the following equivalent integral equations: 


f {t,-l)5uA-[o,.Ju,dQ = Q 


(95) 


(u- - )5t.dr = 0 


(96) 




(97) 


All t and u are unknown. The traction (natural) boundary condition T = C on F, and the 
essential boundary condition u = ;T on I'l, arc satisfied a priori, hence, it can be ignored 
temporarily in the following development. 

It can be seen that equations (95) and (96) hold in any sub-domain. According to 
the concept of the MLPG (Atluri and Shen, 2002a), we use the following weak forms on 
a sub-domain and the corresponding boundary Fs to replace equations (95) and (96) 

l_(t,-r,)v,dr-l<T,,,,-,dQ=o ( 98 ) 
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{. («, - ,7,)(,,rfl- = 0 


(99) 


where v and q arc test functions, and the choice of different test function leads to 
different MLPG methods. If we take the test functions as the weight function in the MLS 
approximation, we will obtain the method developed in Zhang and Yao (2001), however, 
the local weak forms used here arc a little diderent from those in Zhang and Yao (2001), 
where the boundary of the subdomain is involved, that actually should not appear. If 
the subdomain does not intersect with the boundary, equation (98) will not include the 
boundary integral. In equations (98) and (99), u and t on the boundary, and v and q are 
interpolated according to equations (45) and (57b). The u and t inside and on L aie 
defined as [Atluri and Granncll (1978), Zhang and Yao (2001)] 


= S''./'/' 


( 100 ) 


where and are the fundamental solutions with the source point at a node /, a' are 
unknown parameters. We only consider the nodes on the boundary. 

As u is expressed by equation (100), the last integral on the left-hand side of 
equation (98) vanishes if one excludes node J from the subdomain Q.s ^1 which the 
singularity occurs. This singularity will be considered when evaluating the boundary 
integrals. Then by substituting equations (45), (57) and (100) into (98) and (99), and 
omitting the vanished terms, we can obtain the linal system of equations. The equations 
can be solved in the same way as that in Subsection 4.1. 

Similar to Subsection 4.1, four test junctions are chosen to formulate four 
different MLPG methods for boundary variational principle: 
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(1) MLPGl; the test functions v, and t/„ over F, are the same as the weight function in 
the MLS approximation. Then it can be found that we oblian the method 
developed in Zhang and Yao (2001). 

(2) MLPG2; the test functions v, and r/„ over Fy is the collocation Dirac’s Delta 
function (collocation method). 

(3) MLPG5: the test functions v, and r/,, over F., is the Heaviside step function 
(constant over each local sub-domain F,,) 

(4) MLPG6; the test functions v’, and r/„ over F, is identical to the trial function 
(Galerkin method). In this case, let v,- =17, and q-, =/). 


5. Conclusion 

The meshless local Pctrov-Galerkin method (MLPG) method is extended to treat 
the boundary integral equations in this paper. Five boundary integral solution methods 
are introduced; direct solution method; displacement boundary-value problem; traction 
boundary-value problem; mixed boundary-value problem; and boundary variational 
principle. Based on the local weak torm of BIF, four different nodal-based local test 
functions arc selected, leading to lour dillercnt MLPG methods for each BIE solution 
method. These methods combine the advantage of the MLPG and the boundary element 
method. 

A very simple method is used to derive the weakly singular traction boundary 
integral equation based trn the integral relationships for displacement gradients [Okada 
and Atluri, J989]. 

Numerical implementation and demonstrations of the advantages of the presently 
proposed MLPG method foi BIE will be discussed by the authors in forthcoming papers. 
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Captions 



Fig.l Interrelationship of meshlcss methods lor BIE. 
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